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1 Motivation 

[Oct 26, 2005] 

Most of the methods discussed in this course: separation of variables, Fourier 
Series, Green's functions (later) can only be applied to linear PDEs. However, the 
method of characteristics can be applied to a form of nonlinear PDE. 

1.1 Traffic flow 

Ref: Myint-U & Debnath §12.6 

Consider the idealized flow of traffic along a one-lane highway. Let p (x, t) be the 
traffic density at (x,t). The total number of cars in x 1 < x < x 2 at time t is 



Assume the number of cars is conserved, i.e. no exits. Then the rate of change of the 
number of cars in x\ < x < x 2 is given by 



where V (x,t) is the velocity of the cars at (x,t). Combining (1) and (2) gives 




(1) 



— — = rate in at x\ — rate out at x 2 
at 



= p{x x ,t) V (x u t) - p(x 2 ,t) V(x 2 ,t) 




(2) 
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and since xi, x 2 are arbitrary, the integrand must be zero at all x, 

We assume, for simplicity, that velocity V depends on density p, via 

V{p) = c(l--?- 



where c = max velocity, p = p max indicates a traffic jam (V = since everyone is 
stopped), p = indicates open road and cars travel at c, the speed limit (yeah right). 
The PDE (3) becomes 

dt +C { 1 p max ) dx ° ^ 
We introduce the following normalized variables 

p 

u = , t = ct 

Pmax 

into the PDE (4) to obtain (dropping tildes), 

ut + (l-2u)u x = (5) 

The PDE (5) is called quasi-linear because it is linear in the derivatives of u. It 
is NOT linear in u (x,t), though, and this will lead to interesting outcomes. 

2 General first-order quasi-linear PDEs 

Ref: Guenther & Lee §2.1, Myint-U & Debnath §12.1, 12.2 
The general form of quasi-linear PDEs is 

where A, B, C are functions of u, x, t. The initial condition u (x, 0) is specified at 
t = 0, 

u(x,0) = f(x) (7) 

We will convert the PDE to a sequence of ODEs, drastically simplifying its solu- 
tion. This general technique is known as the method of characteristics and is useful 
for finding analytic and numerical solutions. To solve the PDE (6), we note that 

(A,B,a)-( Ux ,ut,-l) = 0. (8) 



2 



Recall from vector calculus that the normal to the surface / (x, y,z) =0 is V/. 
To make the analogy here, t replaces y, f (x, t,z) — u (x, t) — z and V/ = (u t , u x , —I). 
Thus, a plot of z — u (x, t) gives the surface / (x, t, z) = 0. The vector (u x , u t , —1) is 
the normal to the solution surface z = u(x,t). From (8), the vector (A, B, C) is the 
tangent to this solution surface. 

The IC u (x, 0) = / (x) is a curve in the u — x plane. For any point on the initial 
curve, we follow the vector (A, B, C) to generate a curve on the solution surface, 
called a characteristic curve of the PDE. Once we find all the characteristic curves, 
we have a complete description of the solution u (x,t). 



2.1 Method of characteristics 

We represent the characteristic curves parametrically, 

x — x (r; s) , t — t (r; s) , u = u (r; s) , 

where s labels where we start on the initial curve (i.e. the initial value of x at t — 0). 
The parameter r tells us how far along the characteristic curve. Thus (x, t, u) are now 
thought of as trajectories parametrized by r and s. The semi-colon indicates that s 
is a parameter to label different characteristic curves, while r governs the evolution 
of the solution along a particular characteristic. 

From the PDE (8), at each point (x, t), a particular tangent vector to the solution 
surface z — u(x,t) is 

(A (x, t,u) ,B (x, t,u) ,C (x, t, u)) . 

Given any curve (x (r; s) ,t(r;s),u (r; s)) parametrized by r (s acts as a label only), 
the tangent vector is 

dx dt <9w N 
dr ' dr ' dr 

For a general curve on the surface z = u(x,t), the tangent vector (A,B,C) will 
be different than the tangent vecto (x r ,t r ,u r ). However, we choose our curves 
(x (r; s) ,t(r;s),u (r; s)) so that they have tangents equal to (A, B, C), 

7T = * ? = ^ = C < 9 > 

or or or 

where (A,B,C) depend on (x,t,u), in general. We have written partial derivatives 

to denote differentiation with respect to r, since x, t, u are functions of both r and 

s. However, since only derivatives in r are present in (9), these equations are ODEs! 

This has greatly simplified our solution method: we have reduced the solution of a 

PDE to solving a sequence of ODEs. 
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Figure 1: Plot of f(x). 

The ODEs (9) in conjunction with some initial conditions specified at r = 0. We 
are free to choose the value of r at t = 0; for simplicity we take r — at t — 0. Thus 
t (0; s) = 0. Since x changes with r, we choose s to denote the initial value of x (r; s) 
along the x-axis (when t — 0) in the space-time domain. Thus the initial values (at 
r = 0) are 

x{0;s) = s, t(0;s) = 0, w(0;s) = /(s). (10) 



3 Example problem 

[Oct 28, 2005] 

Consider the following quasi-linear PDE, 



nil f)n 

^ + (l + CM )^ = 0, «(x,0)=/(x) 



where c = ±1 and the initial condition / (x) is 



1, |x| > 1 
2-|x|, |x| < 1 



1, x < -1 

2 + x, -l<x<0 

2-x, 0<x<l 

1, x > 1 



The function / (x) is sketched in Figure 1. To find the parametric solution, we can 
write the PDE as 

Thus the parametric solution is defined by the ODEs 



dt 
dr 



1, 



dx 
dr 



1 + cu, 



du 
dr 



= 



4 



with initial conditions at r = 0, 



t — 0, x — s, u = u (x, 0) = u (s, 0) = / (s) . 

Integrating the ODEs and imposing the ICs gives 

t (r; s) = r 
u(r;s) = f(s) 

x(r;s) = (l + cf(s))r + s = (l + cf(s))t + s 



3.1 Validity of solution and break-down (shock formation) 

To find the time t s and position x s when and where a shock first forms, we find the 
Jacobian: 



d{r,s) \ t r t s 



dx dt dx dt , „, , . . , „, . . 

Shocks occur (the solution breaks down) where J = 0, i.e. where 

1 

t -- 



cf (s) 

The first shock occurs at 

t s = min 



1 



cf'(s). 

In this course, we will not consider what happens after the shock. You can find more 
about this in §12.9 of Myint-U & Debnath. We now take cases for c = ±1. 

For c = 1, since min/' (s) = — 1, we have 

t 1 - 1 

min/'(s) 

Any of the characteristics where f (s) = min/'(s) = — 1 can be used to find the 
location of the shock at t s = 1. For e.g., with s = 1/2, the location of the shock at 
t s — 1 is 

1\\ 1 ( ( 1\\ 1 



Any other value of s where /' (s) = —1 will give the same x s . 
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For c = —1, since max/' (s) = 1, we have 

t 1 - 1 

s max /' (s) 

Any of the characteristics where /' (s) = max /' (s) = 1 can be used to find the 
location of the shock at t s = 1. For e.g., with s = —1/2, the location of the shock at 
t s — 1 is 

Any other value of s where /' (s) = 1 will give the same x s . 
3.2 Solution Method (plotting u(x,t)) 

Since r = t, we can rewrite the solution as being parametrized by time t and the 
marker s of the initial value of x: 

x(t;s) = (l + cf(s))t + s, u(;s) = f(s) 

We have written u (; s) to make clear that u depends only on the parameter s. In 
other words, u is constant along characteristics! 

To solve for the density u at a fixed time t = to, we (1) choose values for s, (2) 
compute x (t ; s), u(; s) at these s values and (3) plot u (; s) vs. x (t ; s). Since / (s) 
is piecewise linear in s (i.e. composed of lines), x is therefore piecewise linear in s, 
and hence at any given time, u — f (s) is piecewise linear in x. Thus, to find the 
solution, we just need to follow the positions of the intersections of the lines in / (s) 
(labeled by s = —1,0, 1) in time. We then plot the positions of these intersections 
along with their corresponding u value in the u vs. x plane and connect the dots to 
obtain a plot of u (x, t). Note that for c = 1, the s = —1, 0, 1 characteristics are given 
by 

s = -l:a; = (l + c/(-l))*-l = 2f-l 
s = 0:x = (l + cf(0))t + = 3t 
s = l:x = (l + cf(l))t + l = 2t + l 

These are plotted in Figure 2. The following tables are useful as a plotting aid: 



s = 


-1 





1 


u = 


1 


2 


1 


X = 





3 
2 


2 
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t = t s = l 



s = 


-1 





1 


u = 


1 


2 


1 


X = 


1 


3 


3 



A plot of u(x,l/2) is made by plotting the three points (x,u) from the table for 
t = 1/2 and connecting the dots (see middle plot in Figure 3). Similarly, u(x,t s ) = 
u (x, 1) is plotted in the last plot of Figure 3. 

Repeating the above steps for c = — 1, the s = — 1, 0, 1 characteristics are given 

by 

s = -1 : a; = (1 - / (-1)) t - 1 = -1 
s = 0:x = (l-/(0))t + = -f 
s = 1 : x = (1 - /(l))t + l = 1 

These are plotted in Figure 4. We then construct the tables: 



s = 


-1 





1 


u = 


1 


2 


1 


X = 


-1 


i 

2 


1 



t = t s = l 



s = 


-1 





1 


u = 


1 


2 


1 


X = 


-1 


-1 


1 



As before, plots of u (x, 1/2) and u (x, 1) are made by plotting the three points (x, u) 
from the tables and connecting the dots. See middle and bottom plots in Figure 
5. Note that for c = 1 the wave front steepened, while for c = —1 the wave tail 
steepened. This is easy to understand by noting how the speed changes relative to 
the height u of the wave. When c = 1, the local wave speed 1 + u is larger for higher 
parts of the wave. Hence the crest catches up with the trough ahead of it, and the 
shock forms on the front of the wave. When c = — 1, the local wave speed 1 — u is 
larger for higher parts of the wave; hence the tail catches up with the crest, and the 
shock forms on the back of the wave. 

4 Solution to traffic flow problem 

[Oct 31, 2005] 

The traffic flow PDE (5) is 

ut + (l-2u)u x = (11) 
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Figure 2: Plot of characteristics for c = 1. 




Figure 3: Plot of u{x, t ) with c = 1 for t = 0, 0.5 and I. 



and has form (6) with (A, B, C) = (1 — 2u, 1, 0). The characteristic curves satisfy (9) 
and (10) 



dx 
dr 

m 

dr 
du 

dr 



l-2u, x(0) = s, 
1, *(0)=0, 
0, u (0) = f(s). 

Integrating gives the parametric equations 

t = r + ci, u = c 2 , x = (1 - 2u) r + c 3 = (1 - 2c 2 ) r + c 3 
Imposing the ICs gives c± — 0, c 2 = / (s), c 3 = s, so that 

f = r, « = /(s), x = (l-2f(s))r + s = (l-2f(s))t + s (12) 
We can now write 

x(t;s) = (l-2f(s))t + s, u(;s) = f(s) 

Again, the traffic density u is constant along characteristics. Note that this would 
change if, for example, there was a source/sink term in the traffic flow equation (11), 
i.e. 

u t + (1 — 2u) u x = h(x, t, u) 

where h(x,t,u) models the traffic loss / gain to exists and on-ramps at various posi- 
tions. 

4.1 Example : Light traffic heading into heavier traffic 

Consider light traffic heading into heavy traffic, and model the initial density as 

{a, x < 

(l-a)x + a, 0<x<l (13) 
I 1 

where < a < 3/4. The lightness of traffic is parametrized by a. We consider the 
case of light traffic a = 1/6 and moderate traffic a = 1/3. 
From (12), the characteristics are [DRAW] 

(l-2a)t + s, s<0 
x={ (1 - 2a - 2 (3/4- a) s)t + s, 0<s<l 
-t/2 + s, s>l 
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For a — 1/6, we have 



x 



jt + s, s<0 
(l-ls)t + s, 0<s<l 

— |t + S, S>1 



For a — 1/3, we have 



x 



lt + s, s<0 
(|-§s)t + s, 0<s<l 
— |t + s, s>l 

Again, for fixed times t = t , plotting the solution amounts to choosing an appropriate 
range of values for s, in this case — 2 < s < 2 would suffice, and then plotting the 
resulting points u (to, s) versus x (t Q , s) in the rrw-plane. 

The transformation (r, s) — > (x, t) is non-invertible if the determinant of the Ja- 
cobian matrix is zero, 

!M = **( XP ^ )=**( l ~ 2f{s) - 2/,(s)r + 1 U2f( S ),-1^0. 
d(r,s) \t r t s J \ 1 J 

(14) 

Solving for r and noting that t — r gives the time when the determinant becomes 
zero, 

t = r= W¥Y (15) 

Since times in this problem are positive t > 0, then shocks occur if /' (s) > for some 
s. The first such time where shocks occur is 

tshock = 2max{f (s)Y (16) 

In the example above, the time when a shock first occurs is given by substituting 
(13) into (16), 

1 1 

tshock 



2 max {/'(*)} 2 

Thus, lighter traffic (smaller a) leads to shocks sooner! The position of the shock at 
tshock is given by 



\ ~ a 

X shock (1 2cCj tshock "q ■ 

i - at 
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